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INTRODUCTION 


Turbulence ingestion noise is a significant part of the broadband noise produced by heli- 
copters. This sound generation mechanism involves the interaction of the rotating blades 
and the unsteady upwash velocity vectors associated with the atmospheric turbulent ed- 
dies ingested by the rotor. The resulting noise exists in the presence of loading noise and 
thickness noise mechanisms. It is also additive to the rotor trailing edge noise and vortex 
interactions noise generation mechanisms which represent additional broadband acoustic 
sources. 

This module incorporates three separate models: an atmospheric turbulence model, a 
rapid distortion turbulence contraction model, and a non-isotropic turbulent inflow acous- 
tic source model. Together they determine the atmospheric turbulence ingestion noise 
produced by a helicopter in flight. 

Inputs to the combined mean inflow and turbulence models are controlled by atmo- 
spheric wind characteristics and helicopter operating conditions. 
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LIST OF SYMBOLS AND NOMENCLATURE 


A Constant in geostrophic drag law 

B Constant in geostrophic drag law 

Ct Rotor thrust coefficient = T/pirR 2 (ClR) 2 
d Zero plane displacement height, m, (ft) 

£ Unit vector along one of the three upstream or downstream axes 

E(k) Wavenumber energy spectrum = :)dcr 
f c Coriolis parameter, 2 u> sin 9 

g Gravitational acceleration, m/s 2 , (ft/s 2 ) 

G Geostrophic wind speed, m/s, (ft/s) 

H General rooftop level, m, (ft) 

k Wave number, radians/m, (radians/ft) 

k von Karman constant 

k Wavevector of turbulence 

kg Wave number defined by Equation 21 

L Monin-Obukhov length = ~^ =£ 3 , m, (ft) 

Lw X Integral length scale = / 0 °° R ww (r x )dr x , m, (ft) 
n Direction of the principle normal to the streamline 

p Exponent in mean velocity power law 

r x Displacement in x direction, m, (ft) 
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Rotor radius, m, (ft) 

Cross correlation coefficient = 

Richardson number = 

Temperature, K, (R) 

Thrust, N, (lbf) 

Turbulence velocity fluctuation, m/s, (ft/s) 
Local time mean velocity, m/s, (ft/s) 
Horizontal freestream velocity, m/s, (ft/s) 
Friction velocity = yr«,/p, m/s, (ft/s) 

Vertical freestream velocity, m/s, (ft/s) 

Rotor induced velocity = m/s, (ft/s) 

Vertical velocity fluctuations, m/s, (ft/s) 

RMS turbulence intensity 
Upstream Cartesian coordinate system 
Height above ground, m, (ft) 

Roughness length, m, (ft) 

Tropopause height, m, (ft) 

Rotor tip path plane angle of attack, degrees 

Gamma function 

Boundary layer thickness, m, (ft) 


3 



0 


© 


A 

A 

P 


£,■ 


T 


T w 


$ 


m 


*«(*) 


x 

n,u> 


w,n 


Alternating tensor 
Geographic latitude, degrees 
Potential temperature, K, (R) 

Wavelength of Fourier component of turbulence, m, (ft) 

Rotor inflow ratio, = Usa 
Rotor advance ratio = Ua °nR ~ 

Density, kg/m 3 , ( Ibm/ft 3 ) 

Downstream Cartesian coordinate system 
Temperature fluctuations, K, (R) 

Wall shear stress, N/m 2 , ( Ibf / ft 2 ) 

Nondimensional wind speed gradient = 

Spectrum tensor = exp(—tk • r)Ru(r)dr 

Wake skew angle = arctan y 0 ] > degrees 

Rotation rate, radians/s 

Vorticity field and magnitude of a Fourier component of turbulence. 
See Equations 24 and 25 . 


Superscripts 

U and D Specify the upstream or downstream location 
Subscripts 

ij,k Specify either a vector (such as one of the vectors £1,62, £3) or a component 

of a vector (such as the first Cartesian component of the vector £2) 
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ATMOSPHERIC TURBULENCE PROPERTIES 


Input 


Name Description 

G Geostrophic wind speed (m/s) 

L Monin-Obukhov stability length (m) 

(note: for neutral conditions, input 0.0) 

0 Geographic latitude (degrees) 

Zo Roughness height (m) 

Z Height above ground (m) 

Output 

Name Description 

G Geostrophic wind speed (m/s) 

Lw x Correlation length scale (m) 

y/w^/U 0 o Vertical component of the rms turbulence normalized by the free stream velocity 

Method 


The atmospheric boundary layer is a fully three dimensional flow with boundary condi- 
tions which change in both space and time. Heat transfer is important as are the Coriolis 
forces due to the earth’s rotation. The complexity of the flow renders a complete ana- 
lytic description of the flow impossible. For these reasons, the description of atmospheric 
turbulence is still basically empirical. 

The model presented below was chosen after a review of the current literature on 
the topic and is viewed as representative of the present state-of-the-art in atmospheric 
turbulence prediction. Most of the correlations are from Snyder [l] . These were chosen 
since they were the only ones found which included non-adiabatic effects (i.e. a non-neutral 
atmosphere where heat transfer is important). 

In order to render the problem tractable, the following simplifying assumptions are 
made: 

1) The flow is stationary. 

2) The flow is homogeneous in any horizontal plane. 


5 




3) The flow is isotropic. 

For neutral conditions, the boundary layer thickness is assumed to be constant at 600 
m. For unstable conditions the boundary layer varies from 1 to 2 km during the course of 
a day. It is assumed here to be constant at 1500 m. For stable conditions, the boundary 
layer thickness is given by the following two equations: 


6 = 


-L + -L- 

30 L .25 U, 


ll 


Z T \ 


( 1 ) 


Equation (1) must be solved iteratively in conjunction with Equation (2) below to 
determine 6. 

The friction velocity, U t , is described by the “geostrophic drag law”: 


In 



= A + In 



(2) 


G, the geostrophic wind speed, is a horizontal equilibrium wind which blows parallel 
to atmospheric isobars. It represents an exact balance between the horizontal pressure 
gradient force and the horizontal component of the coriolis force. Geostrophic wind scales 
give a good approximation to the actual wind in the free atmosphere. The magnitude of 
G is given by: 


G = 



where /„ is the coriolis parameter = 2 tus'mO and V„P is the horizontal pressure gra- 
dient. 

A and B are described by “constants” which are a function of the stability condition 
of the atmosphere: 

l) Neutral Conditions 


A = 1.7 


B = 4. 7 


( 3 ) 
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2) Stable Conditions 


A = ln{6/L) - .96(6 /L) + 2.5 
B = 1.15 (6/L) + 1.1 


( 4 ) 


3) Unstable Conditions 


A = In (-6/L) + {f e 6/V.) + 1.5 

B = T £ 7JT -H.Hf'l/U.)exp(-M/L) (5) 

JcO IU* 

Within the atmospheric boundary layer, there exists a sublayer, called the surface layer, 
where stresses and fluxes are nearly constant. This layer can be as thick as 150 m in neutral 
and unstable conditions, but can be as thin as 10 - 20 m under stable conditions. 

In the surface layer, the log law is valid and a description of the mean velocity profile is 
given by the following equations for neutral, stable, and unstable conditions respectively: 

1) Neutral Conditions 



where d = 0 for Zq < .2 m 
and d = H — Zo/k for Zq > .2 m 

2) Stable Conditions 



3) Unstable Conditions 


( 6 ) 


( 7 ) 


U 

u. 




+ 1/*M* 


+ 2 arctan 


feH) 


( 8 ) 
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where = [l — 15-f ] 4 

Equations (6), (7), and (8) are undefined if Z 0 is zero. For an actual atmospheric 
boundary layer, the surface roughness is never zero. Typical values of surface roughness 
length vary from .01 cm for sand to 100 cm for a pine forest [l] . 

For neutral conditions, the integral length scale is a function of height only and is given 
by: 


L * = AZ 


(9) 


For stable atmospheric conditions, the following equation from Kaimal [2] is used: 


L w x = .015— for .05 < i2, < .2 

Ri 


( 10 ) 


where Ri = 


Richardson number 

ZJL 

1+5Z / L 


Z/L 


for Z/L> 0 
for Zj L < 0 


( 11 ) 

( 12 ) 


For unstable conditions, no correlation is available which takes into account stability 
parameters. Therefore, Equation 9 for neutral conditions is suggested. Strictly, the models 
for the integral length scale are valid only for the surface layer, typically the lower 10 to 
20 percent of the boundary layer. 

Under the assumption of isotropic turbulence, all longitudinal integral scales are equal 
to each other, as are the lateral integral scales, i.e. 


Lu X — Ly Y = Lw Z — Lu (13) 

and 

L{f Y = Ljj Z = Ly X = Ly Z = L]y X — L\y Y — Lv (14) 

Teunissen [3] has shown since all longitudinal correlations are equal and all lateral 
correlations are equal, and by the use of the continuity equation that: 
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Lfj — 2 Ly 


so that 


L v x = 2 L v x = 2L W X (15) 

The vertical component of the rms turbulence intensity is given for neutral, stable and 
unstable conditions respectively by the following equations: 

l) Neutral Conditions 


Vrf r ln(30/Zo) 

U„ • 5P ln(Z/Z 0 ) 


(16) 


where p = .24 + .096 log 10 Zo -f .O16(logi 0 Zo) 2 

2) Stable Conditions (and Neutral) 


Vw* = 1.25f7. for Z/L > -.3 


(17) 


3) Unstable Conditions 

) U, for Z/L < -1 

Since these correlations are strictly valid only for the surface layer, following Couni- 
han’s [4] suggestion a linear interpolation is used between the value at the top of the surface 
layer (assumed to be 15% of the boundary layer thickness) to a value of .01 at the top of 
the boundary layer. 

For neutral conditions, the longitudinal and lateral turbulence intensities are approxi- 
mately: 




y/v? = 2Vty* 
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and 


= l.sVtU^ 


For nonadiabatic atmospheric boundary layers, the horizontal velocity components do 
not obey Monin-Obukhov similarity (as the vertical velocity fluctuations do) and no simple 
formula to describe them exists at the present. 

The von Karman model is used to describe the energy spectrum for isotropic turbulence: 


E(k) = 



(19) 


where 


55 T(5/6) u 2 
<\/^r(i/3) k E 5 


6.25278u 2 L 5 


(20) 


and 


, _ V^r(5/6) .746834 

E ~ L r(l/3) “ L 


( 21 ) 
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DISTORTION TENSOR CALCULATION 


Input 


Note that either meters or feet can be used for input as long as the user is consistent 
throughout. 


NAME 

NSS 

NISO 

ISTR 


ITIME 


XMIN 

XMAX 

YMIN 

YMAX 

ZMIN 

ZMAX 

V 

SIR 

a 

V 0 

R 

Voo 

RERROR 

AERROR 

dS 


FUNCTION 

Relative time step resolution (s) 

Number of streamlines to compute 

Flag for outputing a streamline for plotting 

0 for no output file 

1 for streamline file 

Flag for outputing a deformation file 

0 for no output file 

1 for deformation tensor file 

(Note: either ISTR or ITIME can be 1 but not both at the same time) 
Minimum streamline domain limit in X direction (program will stop 
if it attempts to compute outside of this domain) 

Maximum streamline domain limit in X direction, radii 
Minimum streamline domain limit in Y direction, radii 
Maximum streamline domain limit in Y direction, radii 
Minimum streamline domain limit in Z direction, radii 
Maximum streamline domain limit in Z direction, radii 
Rotor advance ratio 

Rotor tip speed (rotational frequency time radius) 

Rotor tip path plane angle of attack, degrees 
Induced velocity at the rotor, m/s, (ft/s) 

Rotor radius, m, (ft) 

V component velocity to be added to the hover case to simulate 
vertical ascent or descent, m/s, (ft/s) 

Relative error allowed in ODE solver 
(recommended value = 1.0E-10) 

Absolute error allowed in ODE solver 
(recommended value = 1.0E-10) 

Delta distance between streamlines used to compute the deformation 
tensor, radii 

(recommended value = 1.0E-4) 
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RECPOL 


NX 

NY 

NZ 

XI 

X2 

Yl 

Y2 

Z1 

Z2 

NR 

NT 

R2 

©i 

©2 


Coordinate specification trigger 

’P’ for polar specification of streamline starting coordinates 

‘C‘ for cartesian specification of streamline starting coordinates 

Number of X coordinate specification points 

Number of Y coordinate specification points 

Number of Z coordinate specification points 

Minimum X for cartesian specification, radii 

Maximum X for cartesian specification, radii 

Minimum Y for cartesian specification, radii 

Maximum Y for cartesian specification, radii 

Minimum Z for cartesian specification, radii 

Maximum Z for cartesian specification, radii 

Number of radial specification points 

Number of angular specification points 

Maximum radius for polar specification 

Mimimum angle for polar specification, degrees 

Maximum angle for polar specification, degrees 


Output 


NAME FUNCTION 

TITLE 60 character alphanumeric title 
NR Number of radial specification points 

NT Number of angular specification points 

IH Code for selection of nonhomogeneous case 

T Vortex circulation strength, m 2 /s, ( ft 2 /s ) 

TSAV Streamline drift time, s 

(time for fluid particle to travel from upstream to rotor) 

U Velocity component in x direction, radii/sec 

V Velocity component in y direction, radii/sec 

W Velocity component in z direction, radii/sec 

XTSAV Downstream X coordinate of streamline in rotor plane, radii 
YTSAV Downstream Y coordinate of streamline in rotor plane, radii 
ZTSAV Downstream Z coordinate of streamline in rotor plane, radii 
XSLl Downstream X coordinate of streamline in standard cartesian 
coordinate system, radii 

(standard coordinate system has x,y plane parallel to horizon with x axis in flight direct 
YSL1 Downstream Y coordinate of streamline in standard cartesian 
coordinate system, radii 
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ZL1 

XSL2 

YSL2 

ZSL2 

USAV 
VSAV 
WSAV 
X1A1, ETC 


Downstream Z coordinate of streamline in standard cartesian 
coordinate system, radii 

Upstream X coordinate of streamline in standard cartesian 
coordinate system, radii 

Upstream Y coordinate of streamline in standard cartesian 
coordinate system, radii 

Upstream Z coordinate of streamline in standard cartesian 

coordinate system, radii 

Upstream U velocity component, radii/sec 

Upstream V velocity component, radii/sec 

Upstream W velocity component, radii/sec 

Deformation tensor 


If ISTR = 1, a file is created named SFILE which contains the coordinate specifications 
of the streamlines 


The output parameters are described below. 

TITLE 60 character alphanumeric title 

NISO Number of streamlines computed 

NSS Relative time step resolution 

XMIN Minimum X streamline domain limit, radii 

XMAX Maximum X streamline domain limit, radii 

YMIN Minimum Y streamline domain limit, radii 

YMAX Maximum Y streamline domain limit, radii 

ZMIN Minimum Z streamline domain limit, radii 

XPLOT X streamline coordinate, radii 

YPLOT Y streamline coordinate, radii 

ZPLOT Z streamline coordinate, radii 
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Method 


Mean Flow Model 

To calculate the distortion of the atmospheric turbulence field as it convects towards 
the helicopter rotor, it is first necessary to describe the flow into the rotor. The mean 
flow field is needed to define the streamlines (particle paths), and these in turn are used 
in computing the change in vorticity and wavenumber between far upstream and the rotor 
disk. A necessary part of the mean flow model is thus a procedure for computing these 
particle paths. 

The method chosen to calculate the mean flow into the helicopter rotor is to model 
the helicopter wake by a series of ring vortices. A graphical representation of the model is 
shown in Figure 1. The overall flow is made up of the velocity field due to the superposition 
of twenty ring vortices plus an ambient velocity. The spacing is varied so that the vortices 
are placed closer together near the rotor disk (see Figure 1), with the circulation per 
unit length (along the wake) constant. The vortices extend to approximately ten rotor 
diameters along the line of the wake. 

The wake skew angle, X, (the angle between the normal to the rotor disk and the wake 
boundary) is defined by : 


X = arctan 


— Vcosa 
V sin a. — V 0 . 


( 22 ) 


The strength of the vortices is found by matching the predicted induced velocity at the 
rotor face with the input value. 

From momentum theory for helicopter aerodynamics, the rotor induced velocity is: 


r _ URCt 
0 = 2vV ! + A ! 


(23) 


where 


H = Rotor advance ratio = 


{/<*, cos a 
UR 


and 


A = Rotor inflow ratio = 


Uoo sin a — V 0 
UR 
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As stated, the rotor wake is described by a series of twenty ring vortices, along the 
skewed centerline of the wake. The velocity field at any point is found from the summation 
of the velocity field due to the individual vortices. The potential velocity for a ring vortex 
is given by Kuchemann and Weber [5] and by Castles and De Leeuw [6]. The axial and 
radial velocity at any point may be expressed by: 

v - = ^ AB+CD v 


V R = 


-r 

2wxR 


{AB + CDF) 


where 


A = K{t) - E(r) 


B = 


x - 1 x + 1 
d\ d,2 


C — di -f* 


D = 


tE(t) 
1 — r 2 


F = 


1 + s 2 + z 2 - djdj 

2x 2 


(l + x)d t 2 — (1 — x)d 2 2 
2xdid,2 



' z 
F - - 
x 


l + x 2 + z 2 

d\d,2 
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di = \fz 2 + (x - l ) 2 


d 2 = \J z 2 + [x + l ) 2 


0^2 — ^1 
r = — 

®2 " I " ®1 

x = radial distance from the axis of the vortex ring nondimensionalized by the ring 
radius 

z = axial distance from the plane of the vortex ring nondimensionalized by the ring 
radius 


K(t) = complete eliptic integral of the first kind 

rx/2 1 


= /' 


\/l — r 2 sin 2 o; 


da 


E(r) 


complete eliptic integral of the second kind 



Given the helicopter operating parameters, the procedure described above can be used 
to compute the velocity at any station in the flow field. A streamline can then be computed 
by integrating the streamline equations: 


dx ■ = Udt 
dy = V dt 

dz = Wdt 
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where U, V, and W are the local time mean velocities in the x, y, and z directions. 

These are three, coupled, ordinary differential equations. Starting at a given station, 
these expressions can be integrated forward or backward in the time domain to predict the 
next streamline position. A computer subroutine, written by Shampine and Gordon [7] 
was used for numerical integration. This method uses a modified divided difference form 
of the Adams - Pece formulas and local extrapolation. It adjusts the order and step size 
to control the local error per unit step in a generalized sense. 


Turbulence Contraction Model 

A “Rapid Distortion Theory” approach has been adopted for calculating the evolution 
of the turbulent flow. Basically the analysis models the flow field in terms of a primary 
flow, induced by the rotor, which convects (and distorts) the vortex filaments which are 
used to represent the turbulence. The rapid distortion approach thus accounts for the 
inherently three-dimensional processes of vortex stretching and tilting, which can change 
both the frequency distribution and the intensity of the turbulence. 

The assumptions inherent in this approach are: 1) the turbulence intensity is small, 2) 
the overturning time for a turbulent eddy is much larger than the time for the distortion 
of the mean flow, and 3) the Reynolds’ number associated with the turbulence is large. 

The basic concept of the Rapid Distortion Theory is to relate the components of the 
vorticity at given locations along a streamline. In a linearized rapid distortion analysis, 
these streamlines (along which the fluid particles are tracked) are streamlines associated 
with the mean flow. Thus the vortex filaments are regarded as being convected by the 
mean flow. The velocities induced by the perturbation vorticity are therefore neglected in 
considering the deformation of the vortex filaments. 

A general turbulent velocity field for incompressible flow can be expressed as a distri- 
bution of vorticity superimposed on a potential flow field. For an inviscid constant density 
flow, vortex lines move with the fluid in which they are embedded. Thus any deformation 
of the fluid will distort the vortex lines and in a contracting stream, the vortex filaments 
will be stretched and tilted, affecting the velocities induced by the vorticity. 

The distortion of the vortex lines from one point on a streamline to another can be 
described by a deformation tensor. This deformation can be decomposed into a fluid 
rotation (3 quantities), the directions of the principal axes (3 quantities), and the amount of 
contraction or extension along each axis (3 quantities). Denoting the upstream coordinates 
by Xi and the downstream coordinates by , the deformation tensor is denoted by dxi/d£j. 

The basic equation for the vorticity at any location in terms of the initial vorticity 
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contained by a given fluid particle is: 


( 24 ) 


( 25 ) 

The deformation tensor is calculated using the streamline tracing procedure described 
above. Starting at a given point on the rotor face the equations describing the streamline 
coordinates are numerically integrated to trace a streamline from the rotor face to a far 
upstream location, usually taken as three radii from the rotor hub. A second streamline 
starting a distance d£i from the first is then also traced for the same length of time as the 
first. At the upstream end of these streamlines a vector is drawn from the end of the first 
streamline to the end of the second. This vector determines three of the nine components 
of the deformation tensor: dxi/d£i, dx 2 /d£ii dx 3 /d£i. The process is repeated for two 
other streamlines displaced a distance d£ 2 and d£ 3 respectively from the initial one, to 
yield the other six components of the distortion tensor. 


or 


n.v _ n POli 
», -II, 


n, D = 

ax. 


i 


I 
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